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Abstract. In this paper we apply the recently developed mimetic discretization method to 
the mixed formulation of the Stokes problem in terms of vorticity, velocity and pressure. The 
mimetic discretization presented in this paper and in 50 is a higher-order method for curvilinear 
quadrilaterals and hcxahcdrals. Fundamental is the underlying structure of oriented geometric 
objects, the relation between these objects through the boundary operator and how this defines 
the exterior derivative, representing the grad, curl and div, through the generalized Stokes 
theorem. The mimetic method presented here uses the language of differential fc-forms with k- 
cochains as their discrete counterpart, and the relations between them in terms of the mimetic 
operators: reduction, reconstruction and projection. The reconstruction consists of the recently 
developed mimetic spectral interpolation functions. The most important result of the mimetic 
framework is the commutation between differentiation at the continuous level with that on the 
finite dimensional and discrete level. As a result operators like gradient, curl and divergence are 
discretized exactly. For Stokes flow, this implies a pointwise divergence-free solution. This is 
confirmed using a set of test cases on both Cartesian and curvilinear meshes. It will be shown 
that the method converges optimally for all admissible boundary conditions. 



1. Introduction 

We consider Stokes flow, which models a viscous, incompressible fluid flow in which the inertial 
forces are negligible with respect to the viscous forces, i.e. when the Reynolds number is very small, 
Re <C 1. Since Re — UL/v, small Reynolds numbers appear when either considering extremely 
small length scales, when dealing with a very viscous liquid or when one treats slow flows. Despite 
the simple appearance of Stokes flow model, there exists a large number of numerical methods to 
simulate Stokes flow. They all reduce to two classes of either circumventing the Ladyshenskaya- 
Babuska-Brezzi (LBB) stability condition or satisfying this condition, [35]. The first class can 
roughly be split into two subclasses, one is the group of stabilized methods, see e.g. QTJ [55] and 
the references therein, the other the group of least-squares methods, see e.g. 13, 46 . 

The class that satisfies the LBB condition is the group of compatible methods. In compatible 
methods discrete vector spaces are constructed such that they satisfy the discrete LBB condition. 
Best known are the curl conforming Nedelec [55 and divergence conforming Raviart-Thomas [64] 
and Brezzi-Douglas-Marini |21| spaces. A subclass of compatible methods consists of mimetic 
methods. Mimetic methods do not solely search for appropriate vector spaces, but try to mimic 
structures and symmetries of the continuous problem, see [TH [501 El \EH HOI HJ • As a con- 
sequence of this mimicking, mimetic methods automatically preserve structures of the continuous 
formulation. 

At the heart of the mimetic method we present is the generalized Stokes theorem, which cou- 
ples the exterior derivative to the boundary operator. In vector calculus this theorem is equivalent 
to the classical Newton-Leibnitz, Stokes circulation and Gauss divergence theorems. These well- 
known theorems relate the vector operators grad, curl and div to the restriction to the boundary of 
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a manifold. Therefore, obeying geometry and orientation will result in satisfying exactly the men- 
tioned theorems, and consequently performing the vector operators exactly in a finite dimensional 
setting. This is indeed what we are looking for and what our mimetic method has in common with 
finite volume methods, [22 EH ■ In a three dimensional space we distinguish between four types 
of submanifolds, that is, points, lines, surfaces and volumes, and two types of orientation, namely, 
outer- and inner-orientation. The inner and outer orientations can be seen as generalizations of 
the concept of tangential and normal in vector calculus, respectively. This geometric structure 
will form the backbone of the mimetic method to be discussed in this paper. It will reappear 
throughout the paper in various guises. Examples of submanifolds in R 3 are shown in Figure [l] 
together with the action of the boundary operator. 



Outer Orientation 




FIGURE 1. The four geometric objects possible in M 3 , point, line, surface and vol- 
ume, with outer- (above) and inner- (below) orientation. The boundary operator, 
d, maps fc-dimensional objects to (k — l)-dimensional objects. 

By creating a quadrilateral or hexadedral mesh, we divide the physical domain in a large number 
of these geometric objects, and to each geometric object we associate a discrete unknown. This 
implies that these discrete unknowns are integral quantities. Since the generalized Stokes theorem 
is an integral equation, it follows for example that taking a divergence in a volume is equivalent 
to taking the sum of the integral quantities associated to the surrounding surface elements, i.e. 
the fluxes. So using integral quantities as degrees of freedom to perform a vector operation like 
grad, curl or div, is equivalent to taking the sum of the degrees of freedom located at its boundary. 
These relations are of purely topological nature and can be seen as the horizontal connections 
between the geometric objects in Figure [T] The vertical connections - not shown in Figure [T]-, 
however describe the metric dependent parts, which are better known as the constitutive relations. 

In this work we use the language of differential geometry to identify these structures, since it 
clearly identifies the metric and metric-free part of the PDEs. The latter has a discrete counter- 
part in the language of algebraic topology. In mimetic methods we employ commuting diagrams 
to indicate the strong analogy between differential geometry and algebraic topology. The most 
important commuting property employed in this work is the commutation between the projection 
operator and differentiation in terms of the exterior derivative. This means that also in finite 
dimensional spaces, operations like gradient, curl and divergence are performed exactly. This im- 
plies, among others, and most importantly that incompressible Navier-Stokes and Stokes flow are 
guaranteed to be pointwise divergence-free, because the projection operator commutes with the 
divergence operator. 

The similarities between differential geometry and algebraic topology in physical theories were 
first described by Tonti, [71] . A mimetic framework relating differential forms and cochains was 
initiated by Hyman and Scovel, [45], and extended first by Bochev and Hyman, [14], and later 
by Kreeft, Palha and Gerritsma [50]. A framework, closely related to the mentioned mimetic 
framework, is the finite element exterior calculus framework by Arnold, Falk and Winther [5J [5] . 
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A more geometric approach is described in the work by Desbrun et al. [551 US]- An excellent 
introduction and motivation for the use of differential forms in the description of physics and the 
use in numerical modeling can be found in the 'Japanese papers' by Bossavit, [T51 ITS] , 

We make use of spectral element interpolation functions as basis functions. In the past nodal 
spectral elements were mostly used in combination with Galerkin (GSEM) p~Ql |47], and least- 
squares formulations (LSSEM) 59, [HI]- The GSEM satisfies the LBB compatibility condition by 
lowering the polynomial degree of the pressure by two with respect to the velocity. This results 
in a method that is only weakly divergence-free, meaning that the divergence of the velocity field 
only convergence to zero with mesh refinement. The LSSEM circumvents the LBB condition in 
order to be able to use equal order polynomials. The drawback of this method is the poor mass 
conservation property, (3H1 E2] ■ 

The present study uses mimetic spectral element interpolation or basis functions on curvilinear 
quadrilaterals and hexahedrals of arbitrary order as described in J37J HO]- The mixed mimetic 
spectral element method (MMSEM) satisfies the LBB condition and gives a pointwise divergence- 
free solution for all mesh sizes. The mimetic spectral element interpolation functions are tensor 
product based interpolants. In every coordinate direction either a nodal or an edge interpolation 
function is used. By using tensor products, we are able to interpolate points, lines, surfaces, 
volumes, hyper-volumes and higher degree n-cube manifolds. 

Although mimetic spectral elements are used to simulate Stokes flow and to derive numeri- 
cal properties, alternative compatible/mimetic functions could be used in combination with the 
mimetic framework without much change, e.g., compatible B-splines, |23 [ 1241 l33j. and mimetic 
B-splines, [7J. 

This paper is organized as follows: first in Section [2] the Stokes problem in terms of vector 
calculus is given, with its relation to geometry and orientation. In Section[3]a brief summary of the 
most important concepts from differential geometry is given. Section [3] discusses the discretization 
of the Stokes model. It introduces the discrete structures of algebraic topology and a set of mimetic 
operators relating differential forms to cochains; the reduction operator, 1Z, the reconstruction 
operator, X, and its composition, the bounded cochain projection, it^ :=XoTZ. As reconstruction 
functions the mimetic spectral element basis functions arc used in this paper. A mixed formulation 
for the Stokes problem is formulated in Section [5] In Section [6] numerical results are discussed 
that show optimal convergence of all variables on curvilinear quadrilateral meshes. Secondly, the 
lid-driven cavity problem is shown on a square, cubic and triangle domain. The last test case is 
the flow around a cylinder moving with a constant velocity. 

2. Stokes problem in vector calculus 

Let n C R™, n — 2, 3, be a bounded n-dimensional domain with boundary dQ. On this domain 
we consider the Stokes problem, consisting of a momentum equation and the incompressibility 
constraint, resulting from the conservation of mass. The Stokes problem is given by 

(2.1a) V • a = f on ft, 

(2.1b) divw = on Q, 

where the stress tensor a is given by 

(2.2) a = -uVu+pI, 

with u the velocity vector, p the pressure, / the forcing term and v the kinematic viscosity. In 
case of velocity boundary conditions the pressure is only determined up to a constant. So in a 
post processing step either the pressure in a point in can be set, or a zero average pressure can 
be imposed; i.e. 

(2.3) / P dn = o. 

Ja 

For the method we like to present, we want to restrict ourselves to vector operations only. 
Therefore, instead of considering the divergence of a stress tensor, V • (z/Vu), we write this as 
uAu by considering constant viscosity. Then the following vector identity is used for the vector 
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Laplacian, — Au = curl curl u — grad div u. The vorticity- velocity-pressure formulation is obtained 
by introducing vorticity as auxiliary variable, Co — curlu. In terms of a system of first-order PDEs, 
the Stokes problem becomes 

(2.4a) to — curl it = on fi, 

(2.4b) curlw + gradp = / on fl, 

(2.4c) divu = onfl. 

Since these PDEs should hold on a certain physical domain, we will include geometry by means 
of integration. In that case we can relate every physical quantity to a geometric object. Starting 



with the incomprcssibility constraint (2.4c) we have due to Gauss' divergence theorem, 



divitdV = / u-ridS = Q, 
v Jdv 



and by means of Stokes' circulation theorem the relation (2.4a) can be written as 



uj x ridS = / curlwxnd<S= / u-tdl. 

S JS JdS 

From the first relation it follows that div it is associated to volumes. The association to a geometric 
object for velocity u is less clear. In fact it can be associated to two different types of geometric 
objects. In the incompressibility constraint velocity denotes a flux through a surface that bounds a 
volume, while in the circulation relation velocity is defined along a line that bounds the surface. We 
will call the velocity vector through a surface outer- oriented and the velocity along a line segment 
inner- oriented. Similarly, vorticity has two different representations, either as the rotation in a 
plane as Stokes' circulation theorem above suggests, or the Biot-Savart description of rotation 
around a line. In the former case Co is inner- oriented whereas in the latter case Co is outer- oriented, 
see also Figure [Tj In fact both the velocity vector u and the vorticity vector Co itself do not have a 
connection with geometry. Therefore, it are the terms u-tdl, u ■ ndS, Co x ndS and Co x t dl that 
are more useful variables when considering Stokes problem on a physical domain. 



The last equation to be considered is (2.4b). This equation shows that classical Newton- 



Leibnitz, Stokes circulation and Gauss divergence theorems tell only half the story. From the 
perspective of the classical Newton-Leibnitz theorem, the gradient acting on the pressure relates 
line values to their corresponding end point, while the Stokes circulation theorem shows that the 
curl acting on the vorticity vector relates surface values to the line segment enclosing it. So how 
does this fit into one equation? In fact there exists two gradients, two curls and two divergence 
operators. One of each related to the mentioned theorems as explained above. The others are 
formal adjoint operators, so the second gradient is the adjoint of a divergence that is related to 
Gauss divergence theorem, the second curl is the adjoint of the curl related to Stokes circulation 
theorem and the second divergence is the adjoint of the gradient related to the classical Newton- 
Leibnitz theorem. Let grad, curl and div be the original differential operators associated to the 
mentioned theorems, then the formal Hilbert adjoint operators grad*, curl* and div* are defined 
as, 

(a, grad* 6) := (diva, 6), (a, curl* 5) := (curia, 6), (a, div* 6) :— (grad a, b). 

Adjoint operators relate geometric operators in opposite direction. Where div relates a vector 
quantity associated to surfaces to a scalar quantity associated to a volume enclosed by these 
surfaces. Its adjoint operator, grad*, relates a scalar quantity associated with a volume to a 
vector quantity associated with its surrounding surfaces. This is illustrated in Figure [2] Following 
Figure [2| the adjoint operators consist of three consecutive steps: First, switch to the other type 
of orientation (inner — > outer or outer —¥ inner), then take the derivative and finally switch the 
result back to its original orientation. 



So (2.4b) could then either be associated to an inner-oriented line segment by rewriting it as 

curl* Co + gradp = /, 
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Figure 2. Geometric interpretation of the action of the boundary operators, 
vector differential operators and their formal Hilbert adjoint operators. 



or be associated to an outer-oriented surface by rewriting it as 

curia; + grad*p = /. 

Without geometric considerations we could never make a distinction between grad, curl and div 
and their associated Hilbert adjoints div*, curl* and grad*. Vector calculus does not make this 
distinction. 

Since our focus is on obtaining a pointwise divergence-free result, we decide to use a formulation 
associated to outer-oriented geometric objects. Then the Stokes problem becomes, 

(2.5a) lu — curl* u = 0, 

(2.5b) curl oj + grad* p = f, 

(2.5c) divu = 0, 

where the first equation is associated to line segments, the second to surfaces and the third to 
volumes. In |12l 113] the same velocity-vorticity-pressure formulation is given in terms of grad, 
curl, div and grad*, curl* and div*. 

For a valid equation, the mathematical objects should be the same; we can only add vectors 
with vectors and scalars with scalars, but not scalars with vectors. But now we add that equations 
also need the be geometrically compatible. We can only add quantities associated with the same 
kind of geometry and with the same type of orientation. This lack of geometric notion in vector 
calculus is what motivates many to use the language of differential geometry instead, H3 [3 
HH [TH [HI EH EH EH SSI SSI [SOI IZl]. Other advantages of using differential geometry over vector 
calculus are that it possesses a clear distinction between variables associated with inner- and outer- 
orientation and it makes a clear distinction between topological and metric-dependent operations. 
All horizontal realtions in Figure [2] are topological. Any detour along geometric objects with 
the other type of orientation introduces metric in the equation. In differential geometry these 
structures are intrinsically embedded. It naturally leads to a discretization technique that can be 
seen as a hybrid between the finite volumes (topological part) and finite elements (metric part). 

3. Differential geometry 

This section presents the Stokes model in the language of differential forms. Differential geom- 
etry offers significant benefits in the construction of structure-preserving spatial discretizations. 
For example, the generalization of differentiation in terms of the exterior derivative encodes the 
gradient, curl and divergence operators from vector calculus and the codifferential represents the 
associated Hilbert adjoint operators grad*, curl* and div*. The generalized Stokes theorem encap- 
sulates their corresponding integration theorems, respectively. The coordinate-free action of the 
exterior derivative and generalized Stokes theorem give rise to commuting properties with respect 
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to mappings between different manifolds. These kind of commuting properties are essential for 
the structure preserving behavior of the mimetic method. 

Only those concepts from differential geometry which play a role in the remainder of this paper 
will be explained. Much more can be found in [34j l36l [50] . 

3.1. Differential forms. Let A fe (il) denote a space of differential k-forms or k-forms, on a suffi- 
ciently smooth bounded n-dimensional oriented manifold il C K™, x := (or, . . . ,x n ), with bound- 
ary d£l. Every element € A fc (f2) has a unique representation of the form, 

(3.1) a {k) = //(x)da; 11 A da;' 2 A • • • A dx lk , 

i 

where / = i\, . . . and 1 < i\ < . . . < if. < n and where //(x) are continuously differentiable 
scalar functions. Differential forms can be seen as quantities that live under the integral sign, [34 , 
which were indicated in the previous section as the 'more useful variables'. For ft C M 3 with a 
Cartesian coordinate system x :— (x,y,z), the outer-oriented vorticity, velocity and pressure are 

UJ^ = ^i(x) da; + ^(x) dy + W3(x) dz, 
vP^ = u(x) dyAdz + u(x) dzAda; + w(x) dxAdy, 
— p(x) dxAdyAdz. 

We see that ui is associated with line elements da;, dy and dz. This is the outer-oriented represen- 
tation in terms of Biot-Savart of rotation around a line segment. Similarly, velocity is associated 
with surface elements, dyAdz, dzAda;, da:Ady, which is the outer-oriented representation of the 
velocity flux through a surface. Finally, writing pressure as a volume form also corresponds to an 
outer-oriented representation. 

Differential k- forms are naturally integrated over fc-dimensional manifolds, i.e. for 6 A k (Q) 
and fl k C R n , with k = dim(ft fc ), 

(3.2) f a (fc) eM <S> (a (fe) ,n fe ) G K, 

Jn k 

where (•,•) indicates a duality pairing between the differential form and the geometry. This 
duality pairing is a metric- free operation, see |36j . Note that the n-dimensional computational 
domain is indicated as f2, so without subscript. We would like to distinguish between fc-forms that 
can be integrated over outer-oriented fc-dimensional manifolds and fc-forms that can be integrated 
over inner-oriented fc-dimensional manifolds. To emphasis this difference, we sometimes write the 
space of the latter as A k (Q). 

The wedge product, A, of two differential forms and is a mapping: A : A k (il) x A'(f2) — > 
A fc+ '(0), k + l < n. The wedge product is a skew-symmetric operator, i.e. Ab^ = (— l) fe '&(') A 
(jW. The pointwise inner-product of fc-forms, (•, •) : A k (Q) x A k (tt) — > K, is constructed using 
inner products of one- forms, that is based on the inner product on vector spaces, see (341 136j . 

The wedge product and inner product induce the Hodge-* operator, ★ : A fc (£l) — > A n ~ k (n), a 
metric operator that includes orientation, defined as 

(3.3) o<*) A*&W := (a<*>,&< fc V n) > 

where e A™ (ft) is a unit volume form, = *1. Let (da;, dy, dz) be a basis in M 3 for 1-forms 
associated with inner-oriented line segments. Then by applying the Hodge-* we retrieve a basis 
for 2-forms associated with outer-oriented surfaces, 

*da; = dyAdz, -kdy — dzAda;, *dz — dxAdy. 

Therefore, the Hodge operator switches between inner- and outer-orientation. The Hodge-* op- 
eration can be interpreted as the vertical relations as given in Figure [2] and coincides with a 



(•, ■) denotes metric dependent inner products, while (■, ■} denotes metric-free duality pairing. 
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rcf) 



constitutive relation. The space of fc-forms on fl can be equipped with an L 2 inner product, 
(•, -) n : A k (n) x A k (Q) -> R, given by, 

(3.4) {a (k \b^) n := / (a< W° = f a« A 

Jn Jn 

The differential forms live on manifolds and transform under the action of mappings. Let $ : 
fi rc f — > Q be a mapping between two manifolds. Then we can define the pullback operator, 

: A fc (£l) — > A fe (r2 rc f), expressing the fc-form on the reference manifold, ri ro f. The mapping, 
and the pullback, $*, are related by 

(3.5) f a (fc) = f $*a (fc) ^ (a^ k \ $(fi rof )) = ($*a (fc) , ft 

A special case of the pullback operator is the trace operator. The trace of fc-forms to the boundary, 
tr : A fc (ft) —> A fc (<9ft), is the pullback of the inclusion of the boundary of a manifold, dft ft, 
see [50]. 

An important operator in differential geometry is the exterior derivative, d : A fc (ft) — > A fc+1 (ft). 
It represents the grad, curl and div (also rot in 2D) operators from vector calculus. It is induced by 
the generalized Stokes' theorem, combining the classical Newton-Leibnitz, Stokes circulation and 
Gauss divergence theorems. Let ftfc+i be a (fc + l)-dimensional submanifold and e A fc (ft), 
then 

(3.6) f daM=[ tra<*> «■ (da<*>, ft fe+1 ) = <tr a«, dSl k+1 ), 

where dft^+i is a fc-dimensional manifold being the boundary of ftfc+i. Due to the duality pairing 



in (3.6|, the exterior derivative is the formal adjoint of the boundary operator d : ftfc+i — > ftfc as 
indicated by the duality pairing, ( |3.2[ ). The boundary operator defines the exterior derivative. 
The exterior derivative is independent of any metric or coordinate system. Applying the exterior 
derivative twice always leads to the null (k + 2)-form, d(da^ fe ^) = 0^ k+2 K On contractible domains 
the exterior derivative gives rise to an exact sequence, called de Rham complex |36) . and indicated 
by(A,d), 

(3.7) R ^ A°(n) A A x (ft) A ^■A n (O)-^-0. 

In vector calculus a similar sequence exists, where, from left to right for R 3 , the d's denote the 
vector operators grad, curl and div. Both inner- and outer-oriented spaces of differential forms, 
A fc (ft) and A fc (ft), possess a de Rham sequence. The two are connected by the Hodge-* operator, 
and constitute a double de Rham complex, 

R — > A (ft) A A 1 (Q) A ... A A n (ft) A 

(3.8) *X *X *X 

OA A n (ft) A A"- 1 ^) ^~ ■•■ ^~ A°(fl) R. 
Observe the similarity between diagram ( |3.8[ ) and Figures [T] and [2j which is due to the fact that 
the exterior derivative is the adjoint of the boundary operator. The pullback operator and exterior 
derivative possess the following commuting propertjj^J 

(3.9) $*da (fc) = d$*a (fc) , Va (fe) e A fc (ft), 
as illustrated in the following commuting diagram, 

A k (n) — A fe+1 (ft) 

A fc (ft rcf ) — A^ft^). 



2 Note that on the lefthandside of this equation we consider the pullback of a (fc + l)-form, whereas on the right 
hand side the pullback of a fc-form. We could have written this as , j^fcdW = d^^a^). In order to improve 
readibility and knowing that the meaning of these operators is clear from the context we do not explicitely denote 
this. 
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The inner product, (3.4), gives rise to the formal Hilbert adjoint of the exterior derivative, the 
codifferential operator, d* : A fc (^) A*" 1 ^), a s (da**" 1 ), b^) Q = (a^A d*b^) Q , which 
represents the grad*, curl* and div* operators. Whereas the exterior derivative is a metric-free 
operator, the codifferential operator is metric-dependent, and given by d* = (— i)™( fc + 1 )+ 1 * d*, 
[361 150] . Here we see the three operations that were mentioned in the previous section and were 
illustrated in Figure [2j Switch to the other type of orientation, *, apply the derivative, d, and 



switch back to the original orientation, *. In case of non-zero trace, and by combining (3.4) and 



(3.6), we get 



(3.10) (a^ k - 1 \d*b^) n =(da^- 1 \b^) Q - trfAAtr 

Jan 

This is better known as integration by parts and is often used in finite clement methods to avoid the 
codifferential. Also for the codifferential, on contractible manifolds there exists an exact sequence, 

(3.11) A A°(n) A- A x (fi) ^ A A"(fi) ^ R. 

Finally, the Hodge-Laplace operator, A : A k (fl) — > A fc (£l), is constructed as a composition of the 
exterior derivative and the codifferential operator, 

(3.12) - Aa (fe) := (d*d + dd*)a (fe) . 

3.2. Hilbert spaces. Function spaces play an important role in the analysis of numerical meth- 
ods. Of importance in this paper are the Hilbert spaces. On an oriented Riemannian manifold, 



we can defin e Hi lbert spaces for differential forms. Let all //(x) in (3.1) be functions in L (f2) 



then a( fe ) in (3.1 1 is a /c-form in the Hilbert space L 2 A k (fl). The norm corresponding to the space 
L 2 A k (Vt) is ||a^||L2 A k = \J (a^ k \a^)a. Although extension to higher Sobolev spaces are possi- 
ble, we focus here on the Hilbert space corresponding to the exterior derivative. The Hilbert space 
HA k (n) is defined by 

(3.13) HA k (n) = {a<- k) G L 2 A k {n) | da^ e £ 2 A fe+1 (ft)}, 

and the norm corresponding to HA k (fl) is defined as 

(3-14) :-||a( fe )||| 2At + ||d fl ( fe )||i 2Afc+1 . 

The Hilbert complex, (HA, d), a special version of the de Rham complex, is the exact sequence of 
maps and spaces given by 

(3.15) R =-> HA Q {n) A HA^il) A U HA n (Q) A 0. 

In vector operations the Hilbert complex becomes for C R 3 , 

(3.16) H\n) A H(cvsl, n) A 2T(div, fi) A L 2 {ty, 
and for ft C R 2 , either 

(3.17) i?A) Aitf(rot,tt) ^L 2 (!]), or H^Q) A H(div, Q) A L 2 (A 



The two are related by the Hodge-* operator (3.3), see 



HA°(Q) A HA^fl) A L 2 A 2 {Vt) H^Sl) A H(div,Q) A L 2 (Sl) 

(3.18) *$ *$ *$ & *$ *$ *t 

L 2 A 2 (fl) A HA 1 (SI) A HA°(Q) L 2 (n) A iJ(rot,0) £^ fT^fi). 

A similar double Hilbert complex can be constructed in R 3 . Again note the similarities with these 
double Hilbert complexes and that of (3.8) and geometric structure depicted in Figures [l] and [2j 
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3.3. Stokes problem in terms of differential forms. The kind of form a variable has is di- 
rectly related to the kind of manifold this variable can be integrated over. For example, from 
a physics point of view velocity is naturally integrated along a line (streamline), a 1-manifold, 
indicating that velocity is a 1-form. However, looking at the incompressibility constraint, velocity 
in incompressible (Navier)-Stokes equations is usually associated to a flux through a surface, indi- 
cating that velocity should be an (n — l) -form (n = dim(O)). The two are directly related by the 



Hodge duality, ( |3.8[ ). The Hodge-* not only changes the corresponding type 

of integral domain, but also its orientation (along a line = inner, through a surface = outer). 

Note that the Hodge-* is often combined with a constitutive relation. In that case the two 
variables have clearly a different meaning. In incompressible flow models, mass density plays the 
role of material property, so we actually have (pu)^™ -1 ^ = •k p u( 1 '. Since mass density is assumed 
to be equal to one in incompressible (Navier)-Stokes, this difference is less obvious. 

As for the velocity, also for pressure and vorticity there exists an inner and outer oriented 
version. The inner oriented variables are pressure, p G A°(f2), associated to point values, vorticity 
and u) G A 2 (f2), associated to circulation in a surface. Alternatively, there exists the set of 
outer-oriented variables, being the pressure, p G A"(f2), measured in a volume, and vorticity, 
uj G A"~ 2 (f2), corresponding to circulation around a line (both in case of £1 C 

Both sets, (p(°),uW,w( 2 )) and (a/"" 2 ) 
mulations and numerical schemes. For the former see [TJ [THj and for the latter see [HI [31] ■ 

To obtain a pointwise divergence-free solution, the incompressibility constraint is leading, and 
therefore the set of outer-oriented variables are used in this paper, [w^ n ~' i \u^ n ~ 1 ' ,p^ n '), with 
forcing term /t™™ 1 ). Then the Stokes problem in terms of differential forms becomes, 



u ( n 1 ) 5 p(™)) are used in literature to derive mixed for- 



(3.19a) 
(3.19b) 



+ dV n) 



(n-l) 



on f2, 
on O, 



where A is the Hodge-Laplacian defined by (3.12). Vorticity is introduced as auxiliary variable to 



cast this system into a system of first-order equations. Substitution of (3.12) and the incompress 



ibility constraint (3.19b), gives the vorticity- velocity-pressure formulation in terms of differential 
forms, 



(3.20a) 
(3.20b) 
(3.20c) 



Ul 



(n-2) 



i/dw ( ™- 



d * w («-l) 
1 + d*p (n) 

du {r 



0, 

/( „- 



= 0, 



on f2, 
on f2, 
on il. 



Note the resemblance of this system with (2.5). Note also that whereas grad, curl and div are 
only defined in M 3 , (|3.20|) is valid in 1" for all n > 1. 



The actions of the exterior derivatives and codifferentials in this system are illustrated below 
for a two-dimensional domain. 

Example 1 (2D Stokes problem). Let fl C M. 2 , with Cartesian coordinates x := (x,y), and 



let the two-dimensional de Rham complex be equivalent to the second complex in (3.11). Then 
velocity is expressed as 

w 1 ' = — v(x.)dx + u(x)dy. 
Applying the exterior derivative gives us a 2-form, the divergence of velocity, 

du(1) -{^x + W dXAdV - 
Vorticity is a 0-form, u/°) = w(x) G A°(Q), and the curl of vorticity gives, 

d^ = px + py. 
ox ay 



3 With • we indicate a variable contained in the lower complex of (13 . 8b . 
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The gradient of pressure, pfc> — p(x)dx A dy € A 2 (Q), is the action of the codifferential, 



A* (2) d P A , d ^A 

d p y > = -7^-da; + -^-dy. 
ay dx 



Then the momentum equation follows, 



8uj dp\ , / duj <9p\ , ,, , 

dx+ — + — )dy = -f y (x,y)dx + f x (x, y)dy. 



dx dy J \ dy dx 

In a similar way the vorticity-velocity relation can be obtained. 

4. Discretization of Stokes problem 

The mimetic discretization of the Stokes problem consists of three parts. First, the discrete 
structure is described in terms of chains and cochains from algebraic topology, the discrete coun- 
terpart of differential geometry. This discrete structure mimics a lot of properties of differential 
geometry. Secondly, mimetic operators are introduced that relate the continuous formulation in 
terms of differential forms to the discrete representation based on cochains. Thirdly, mimetic 
spectral element basis functions are described which satisfy the structure defined in the algebraic 
topology and mimetic operators sections, Sections |4 . 1 1 and |4~2] respectively. The action of the exte- 
rior derivative, i.e., grad, curl and div, are shown, which leads among others to the divergence- free 
solution. 



4.f . Algebraic Topology. In many numerical methods, especially in finite difference and almost 
all finite element methods, the discrete coefficients are point values, i.e. zero-dimensional sub- 
manifolds. In the structure of algebraic topology, the discrete unknowns represent values on 
fc-dimensional submanifolds, ranging from points to n-dimensional volumes, so < k < n. These 
fc-dimensional submanifolds are called k-cells, r^ k y See [HUSO] [51] how they are formally defined. 
The two most popular classes of fc-cells in literature to describe the topology of a manifold are 
either in terms of simplices, see for instance [54] |68l E3] , or in terms of cubes, see [HI [69] [71] and 
Figure [3] for an example of fc-cubes in K 3 . From a topological point of view both descriptions are 
equivalent, see [30] . Despite this equivalence between simplicial complexes and cubical complexes, 
the reconstruction maps to be discussed in Section [472] differ significantly. For mimetic methods 
based on simplices see [25] UHl HO] ; whereas for mimetic methods based on cubes see [5] HU [SB] ■ 




Figure 3. Example of a 0-cell, a 1-cell, a 2-cell and a 3-cell in K 3 . 



Here we list the terminology to setup a homology theory in terms of n-cubes as given by [51]. 
Consider an oriented unit fc-cube given by I k — I x I x • • • x / (k factors, k > 0), where / = [0, 1] 
is a one-dimensional closed interval. By definition 1° is a space consisting of a single point. Then 
a k-cube in an n-dimensional manifold VI is a continuous map TnA : I k — > tt, < k < n. 

All fc-cells are oriented. This means that we define a default orientation. The default orientation 
of the cell is implied by the orientation of the line segment /, which is defined positive in positive 
coordinate axis direction, and the map T/^y For outer-oriented cells, this for example also implies 
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a positive way of going through a surface and rotating around a line. A fc-cell with opposite 
orientation is said to have a negative orientation. 

The concept of orientation shown in Figures [T] and [2] gives rise to the boundary operator, d, 
that relates a fc-cell to a set of surrounding (fc — l)-cells, which has either the same or opposite 
orientation. Examples are given in Figure [4j where the faces of the fc-cells are shown in black. 




Figure 4. Examples of faces of outer-oriented fc-ceils in M 3 . 

This definition describes the boundary which we already encountered in Figure [T] and ( |3.6[ ). 
The boundary of a fc-cell again consists of a set of (fc — l)-cells, as illustrated in Figure [4] From 
this we can define a cell complex. 

Definition 1 (Cell complex). [42] A cell complex, D, in a compact manifold CI is a finite 
collection of cells such that: 

(1) The set of n- cells in D covers the manifold Q. 

(2) Every face of a cell in D is contained in D. 

(3) The intersection of any two k-cells, Ttfy and ovm in D either share a common I -cell, 
I = 0, . . . , k — 1 in D, is empty, or Tt^ = <T(m . 




O-cells 1 -cells 2-cells 3-cells 



FIGURE 5. Example of a cell complex. Left: a three dimensional compact mani- 
fold. Right: the fc-cells that consistitute the cell complex. 

We call a cell complex an oriented cell complex, once we add to each fc-cell a default orientation 
according to the definition of fc-cubes. Figure [5] depicts an example of a cell complex in a compact 
manifold 51 C 1R 3 . The ordered collection of all fc-cells in D generate a basis for the space of fc- 
chains, Ck{D). Then a fc-chain, £ Ck{D), is a formal linear combination of fc-cells, tqa^ S D, 

(4.1) c {k} =J2°i T (k),i- 

i 
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The fc-cells, TVMj form a basis for the fc-chains. Once such a basis with default orientation has 
been chosen, any chain is completely determined by the coefficients Ci which can be arranged in 
a column vector c = [ci, C2, . . -] T ■ In the description of geometry, we restrict ourselves to chains 
with coefficients in Z/3 = { — 1,0, 1}. The meaning of these coefficients is : 1 if the cell is in the 
chain with the same orientation as its default orientation in the cell complex, -1 if the cell is in 
the chain with the opposite orientation to the default orientation in the cell complex and if the 
cell is not part of the chain. 

We can now extend the boundary operator applied to a fc-cell to the boundary of a fc-chain. 
The boundary operator acting on fc-chains, d : Ck(D) — > Ck-i(D), is defined by [42, 54 , 

(4.2) <9c (fe) = dJ2 c ' T (k)A ■= c%d ( T W,i) ■ 

i i 

The boundary of a fc-cell t^j is a (fc — l)-chain formed by the faces of T(ja. The coefficients of this 
(fc — l)-chain associated to each of the faces is given by the orientations. 

3 

with 

e\ = 1, if the orientation of T(u_i\j equals the default orientation, 

= — 1, if the orientation of Tf k -i),j i s opposite to the default orientation, 
e\ = 0, if T(j._x\j is not a face of T(k),i ■ 

The boundary of a 0-cell is empty. In case all fc-cells in the chain cm have positive orientation, 
so c* = 1, then 

( 4 - 3 ) dc (k) =^2^24nk-i)j- 

i 3 

Recalling that the space of fc-chains is a linear vector space it follows that the boundary operator 
can be represented as a matrix acting on the column vector c of the fc-chain. The coefficients e\ 
are the coefficients of an incidence matrix E( fe -i,fc) that represents the boundary operator. Like 
the exterior derivative, applying the boundary operator twice on a fc-chain gives the null (fc — 2)- 
chain, ddc^ = 0( fe _ 2 ) f° r a U c (k) £ C-k(D), see Figure [61 This was expected, since the exterior 










O 








dc 



/ 



ddc 



FIGURE 6. The boundary of the boundary of a 3-cell is zero because all edges 
have opposite orientation. 



derivative and boundary operator are related according to the generalized Stokes theorem, (3.6). 
This property is reflected in the incidence matrices, since they are matrix representations of the 
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topological boundary operators. Therefore E(k-2.k-i)E(k-i,k) = 0, where for Figure[6]we have 



E(i,: 



-i 
i 
o 
o 
i 

-i 
















1 



-1 



-1 



1 



-(2,3) 



The set of fc-chains and boundary operators gives rise to an exact sequence, the chain complex 
(C k (D),d), 



C k -i{D) *- 



C k (D) f- 



C k+1 (D) +- 



(4.4) . . . <_j 

This sequence is the algebraic equivalent of Figure [l] Dual to the space of fc-chains, Ck(D), is 
the space of k-cochains, C k {D), defined as the set of all linear functionals, c( fe ) : C k {D) 
The duality is expressed using the duality pairing {c^ k \c 



(k)/ 



:=r.( fe )( 



\ k ) 



Note the resemblance 



between this duality pairing and the integration of differential forms, see ( |3.2[ ). 

Let {tvm A } form a basis of Cfc(D), then there is a dual basis {r' fe '' 1 } of C k (D), such that 
T^ k ^ l (r^ i) — Sj and all fc-cochains can be represented as linear combinations of the basis elements, 



(4.5) 



The cochains are the discrete analogue of differential forms. With this duality relation between 
chains and cochains, we can define the formal adjoint of the boundary operator which constitutes 
an exact sequence on the spaces of fc-cochains in the cell complex. This formal adjoint is called 
the coboundary operator, S : C k (D) — > C k+1 (D), and is defined as 

(4.6) (5c« c (fe+1) ) := (c« dc (fc+1) ), Vc« e C k (D) and Vc (fe+1) g C k+1 {D) . 

Also the coboundary operator satisfies SSc^ — o( k+2 ~> for all c( fc ) g C k {D), see Figure [7J and 
gives rise to an exact sequence, called the cochain complex (C k (D),5), 



(4.7) 



-> C k -\D) 



-> C k (D) 



-)• C k+1 {D) 



The coboundary operator is the discrete analogue of the exterior derivative. Also the coboundary 



operator has a matrix representation. As a result of the duality pairing in (4.6), the matrix 



representation of the coboundary operator is the transpose of the incidence matrix of the boun dary 
operator, E( fc ' fc-1) := (E( fc _ l fe) ) T . And again, £( k + 1 , k )£( k < k - 1 ) — o. Note that expression (4.6 1 
is nothing but a discrete generalized Stokes' theorem. The matrices representing the coboundary 
operator only depend on the mesh topology. These matrices will explicitly appear in the final 



matrix system, (5.2 1 




FIGURE 7. The action of twice the coboundary operator S on a 1-cell has a zero 
net result on its surrounding 3-cells, because they all have both a positive and a 
negative contribution from its neighboring 2-cells (reproduced from [14 ). 
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4.2. Mimetic Operators. The discretization of the flow variables involves a projection operator, 
TTh, from the complete space A fc (f2) to a subspace AS (£1; Ck) C A fc (Sl). In this subspace differential 
forms are expressed in terms of fc-cochains defined on fc-chains, and corresponding fc-form interpo- 
lation functions (often called basis- functions) . Usually, the subspace is a polynomial space. The 
projection operation actually consists of two steps, a reduction operator, TZ, that integrates the 
fc-forms on fc-chains to get fc-cochains, and a reconstruction operator, T, to reconstruct fc-forms 
from fc-cochains using the appropriate basis-functions. These mimetic operators were already in- 
troduced before in O 05] • A composition of the two operators gives the projection operator 
-Kh = I o 1Z as is illustrated below. 




K 



C k (D) 

We already saw the similarities between differential geometry and algebraic topology. We now 
impose constraints on the maps 7Z and I to ensure that these structures are preserved. By 
imposing structure-preserving constraints on these operations, these three operators together set 
up the mimetic framework. An extensive discussion on mimetic operators can be found in |50j . 
Here only the most important properties are listed. 

Definition 2 (Reduction). The reduction operator 1Z : A k (tt) —> C k (D) maps differential forms 
to cochains. This map is defined by integration as 



(4.8) 



(IZa 



(fe) 



r(k)) 



a {k) , Vr (fc) e C fc (£»). 



Then for all cr^) € Ck(D), the reduction of the k-form, a^ k ' € A k (Q), to the k-cochain, a.^ £ 
C k (D), is given by 



(4.9) a«(c (fc) ) := (TZa^ 



C (I 



(fc) _ 



7 (fe) 



The reduction map TZ provides the integral quantities that were mentioned in the Introduction. It 
is the integration of a fc-form over all fc-cells in a fc-chain that results in a fc-cochain. A special 
case of reduction is integration of an n-form a € A"(fi) over f2, then 



,(n) ._ 



= (TZ, 



where the chain er 



(n) 



~(ra),i 



(so all c l — +1) covers the entire computational domain f2. The 



reduction map has a commuting property with respect to continuous and discrete differentiation, 
(4.10) Kd = 5U onA k (n). 

This commutation can be illustrated as 

A k _d_^ A fc+i 
C k - 



_> C k+i 

This property follows from the generalized Stokes' theorem (3.6| and the duality pairing of (4.6) 
(4.11) <W fc ),c (fe) ) <U f da« W f a« <U (KaW,dc (k) ) & (61Za< k \ C(fc) >. 

-'c(fc) J dc {k) 

The operator acting in the opposite direction to the reduction operator is the reconstruction 
operator, I. The reconstruction operator I : C k (D) — > A^(Q;Cfc) maps fc-cochains onto finite 
dimensional k- forms. The reconstructed differential forms belong to the space A^il; Ck), which is 
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a proper subset of the complete A:- form space A k (tt). While the reduction step is clearly defined 
in Definition [2j in the choice of interpolation forms there exists some freedom. 

Definition 3 (Reconstruction). Although the choice of a reconstruction method allows for some 
freedom, X must satisfy the following properties: 

• Reconstruction I must be the right inverse oflZ, so it returns identity (consistency prop- 
erty), 

(4.12) TlX = Id on C k (D). 

• Like 1Z, also the reconstruction operator X has to possess a commuting property with respect 
to differentiation. A properly chosen reconstruction operator X must satisfy a commuting 
property with respect to the exterior derivative and coboundary operator, 

(4.13) dX = X6 onC k {D). 
This commutation can be illustrated as 



C k _A_^ c k+1 

Moreover, we want it to be an approximate left inverse of H, so the result is close to identity 
(approximation property) 

(4.14) XK = Id + (h p ) in A fe (ft). 

where 0(h p ) indicates a truncation error in terms of a measure of the mesh size, h, and a polyno- 
mial order p. 

Definition 4 (Projection). The composition X o TZ will denote the projection operator, := 
X1Z : A fe (f2) — » A^(f2;Cfc), allowing for an approximate continuous representation of a k-form 
a W e A k {Vt), 

(4.15) a { h k} = TT h a {k) = XKa {k \ n h a {k) G A k h (n ; C k ) C A k (Q). 

where XTZa^ is expressed as a combination of k-cochains and interpolating k-forms. 

A proof that 7T/j is indeed a projection operator is given in [50] - Since TTha^ = XTZa^ is 
a linear combination of fc-cochains and interpolation fc-forms, the expansion coefficients in the 
spectral element formulation to be discussed in Section 4.3 are the cochains which in turn are the 
integral quantities mentioned in the Introduction. 

Lemma 1 (Commutation property). There exists a commuting property for the projection 
and the exterior derivative, such that 

(4.16) d 7 r fe = 7r/ 1 d on A k (Q). 
This can be illustrated as 



A k 


A k+1 






K - 





Proof. This is a direct consequence of the definitions of the reduction (4.10), reconstruction (4.13) 
and projection operators (4.15), 

dn h a^ & dZ7Za<*> & XSKa^ ® Inda (k) H§ Vo (fc) £ A fe (fi) . 

□ 
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Note that it is the intermediate step X51Za^ k ' that is used in practice for the discretization, 
see Examples [3] and |4j Section |4.4| Since we have a matrix representation of the coboundary 
operator in terms of incidence matrices, we expect the incidence matrices to appear explicitly in 
the spectral element formulation, see (5.2). Lemma [I] is the most important result in this paper. 



As a direct consequence we obtain the pointwise divergence-free solution, as illustrated in the 
following example. 

Example 2. Consider the relation du^™ -1 ) = g^ n ' . In vector notation the d represents the div 
operator. Now let du^ ^ = g^ff 1 be the discretization of our continuous problem. Then by using 



(4.16) we get 

^ ( r 1} - 9 ( h ] = dn h u^ - W n) = 7T h (du^ - gW) = 0. 
It follows that our discretization is exact. In case g^ n ' — 0, we have a pointwise divergence- free 
solution of . 

As a direct consequence of Lemma [l] we satisfy the LBB stability criteria, see [22J, |39j [49] . The 
projection does not commute with codifferential operator. This is the main reason why we rewrite 
the codifferentials into exterior derivatives and boundary integrals, by means of integration by 



parts using (3.10) 



We do not restrict ourselves to affine mappings only, as is required for many other compatible 
finite elements, like Nedelec and Raviart-Thomas elements and their generalizations [SI 1551 Ifi4*]. 
but also allow non-affine maps such as transfinite mappings [40] or isogeometric transformations. 
This allows for better approximations in complex domains with curved boundaries, without the 
need for excessive refinement. This is possible since the projection operator tt^ commutes with 
the pullback $*, 

(4.17) $*Tr h =ir h $* onA fe (n). 

This commutation can be illustrated as 

A k (^) A fc (^rcf) 



A£(fi,C fe ) A*(fi ref ,C fe ) 



An extensive proof is given in 



4.3. Mimetic spectral element basis-functions. Now that a mimetic framework is formulated 
using differential geometry, algebraic topology and the relations between those - the mimetic oper- 
ators - we derive reconstruction functions, I, that satisfy the properties of the mimetic operators. 
In combination with the reduction operator, TZ, it defines the mimetic projection operators, tt/j. 
The finite dimensional A-forms used in this paper are polynomials, based on the idea of spectral 
element methods, Spectral element methods have many desirable features such as arbitrary 
polynomial representation, favourable conditioning, element wise local support, and optimal stabil- 
ity and approximation properties. However, the definition of the reconstruction operator requires 
a new set of spectral element interpolation functions. The mimetic spectral elements were derived 
independently by [371 165] , and are more extensively discussed in [50 . Only the most important 
properties of the mimetic spectral element method are presented here. 

In spectral element methods the computational domain Q is decomposed into M non-overlapping, 
possibly curvilinear quadrilateral or hexahedral, closed sub-domains Q mi 

M 

(4.18) il = \J Q m , Q m n Qi = dQ m n dQ h m ^ I, 



where in each sub-domain a Gauss-Lobatto mesh is constructed, see Figures 10 and 14 in the next 
section. The complete mesh is indicated by Q :— Ylm=i Qm- 

The collection of Gauss-Lobatto meshes in all elements Q m £ Q constitutes the cell complex 
D. For each element Q m there exists a sub cell complex, D m . Note that D m flD|, m ^ I, is not 
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an empty set in case they are neighboring elements, but contains all fc-cells, k < n, of the common 
boundary, see Definition [Tj 

Each sub-domain is mapped from the reference element, Q lc f = [—1,1]", using the mapping 
'I'm : Qref — > Qm ■ Then all flow variables defined on Q m are pulled back onto this reference element 
using the following pullback operation, $* n : A|(Q m ; Ck) — ► A^(Q re f; Cfc). In three dimensions the 
reference element is given by Q re f := {(£, rj, () | — 1 < £, 77, C < 1}- 

The basis-functions that interpolate the cochains on the quadrilateral or hexahcdral elements 
are constructed using tensor products. It is therefore sufficient to derive interpolation functions 
in one dimension and use tensor products afterwards to construct rt-dimensional basis functions. 
A similar approach was taken in [24] . Because the projection operator and the pullback operator 



commute (4.17), the interpolation functions are discussed for the reference element only. 

Consider a 0-form £ A°(Q re f) on Q Te { := £ £ [—1, 1], on which a cell complex D is defined 
that consists of N + 1 nodes, where — 1 < £0 < • • • < £iV < lj an d N edges, rm^ = 
of which the nodes are the boundaries. Corresponding to this set of nodes (0-chains) there exists 
a projection using iV th order Lagrange polynomials, h(£), to approximate a 0-form, as 



N 



(4.19) 7r fc a( > = Y>A(0- 



i=0 



Lagrange polynomials have the property that they interpolate nodal values and are therefore 
suitable to reconstruct the cochain a(°) = IZa^ containing the set ai = a(£j) for i = 0, . . . , N. So 
Lagrange polynomials can be used to reconstruct a 0-form from a 0-cochain. Lagrange polynomials 
are in fact 0-forms themselves, h(£) £ A?(Q re f ; Co). Lagrange polynomials are constructed such 
that their value is one in the corresponding point and zero in all other mesh points, 

(4.20) nim = h(Z P ) = { 1 l]^ p p . 



This satisfies (4.121, where in this case X = k(£). Gerritsma [37] and Robidoux [65] derived a 
similar projection for 1-forms, consisting of 1-cochains and 1-form polynomials, that is called the 
edge polynomial, e^(£) £ A\(Q IC [ ; C\). 

Lemma 2 (Edge polynomial). Following Definitions^ and^ apply the exterior derivative to 
TTha^\ it gives the 1-form ir^b^ = dn^a^ — XSTZa^ given by 

N 

(4.21) ir h b^ =Yj iei (0, 

with 1 -cochain t/ 1 ), where 



=.1 



(4.22) b i = (nb^, T{1) , i )= I 6(0-/ do<°>=/ a' ', 



dT (i),i 



= a (6) - a (6-i) = a i - o-i-1, 
with the edge interpolation polynomial defined by 

i-1 N N i-1 

(4.23) e,-(o = - £ dlk = E d ^ (o = \ E ^ (o - 1 E d ^ (o • 

fc — fc — Z fc — 2 fc — 

Proo/. See [33 [501 [55]. □ 

The value corresponding to line segment (1-cell) mw is given by bi — ai — <2j_i and so b^ = 
(Ja*- ' is the discrete deriva tive operator in ID. This operation is purely topological, no metric is 



involved. It satisfies ( |4.13[ ), since dla (0) = 18sS 0) . Note that we have de^) = J2 do<il i(0 = 0. 
The 1-form edge polynomial can also be written as below, separating the edge function into its 
polynomial and its basis, 

ei(0=£i(0d£, with Ei (£)=-£^. 
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GLL nodal interpolation 




-1 -0.5 0.5 



Figure 8. Lagrange poly- 
nomials on Gauss-Lobatto- 
Legendre mesh. 



GLL edge interpolation 




-1 -0.5 0.5 1 



Figure 9. Edge polynomials 
on Gauss-Lobatto-Legendre 
mesh. 



Similar to (4.20), the edge functions are constructed such that when integrating ej(£) over a line 
segment it gives one for the corresponding element and zero for any other line segment, so 

r£p f 1 if i = p 

(4.24) Keiit) = et{C 



This also satisfies (4.121, where in this case X = The fourth-order Lagrange and third-order 

edge polynomials, corresponding to a Gauss-Lobatto mesh with N = 4, are shown in Figures [S] 
and|U 

Now that we have developed interpolation functions in one dimension, we can extend this to 
the multidimensional framework by means of tensor products. This allows for the interpolation 
of integral quantities defined on fc-dimcnsional cubes. Consider a reference element in IR 3 , Q rc i = 
[— 1, l] 3 . Then the interpolation functions for points, lines, surfaces and volumes are given by, 

point : P$) k &r],0 = k (0 ® h (v) ® h (C) , 

line: 4d,k(^V, = {e 4 (0 ® lj{rj) ® Z fc (C), ® ejfa) ® J fc (C), Z*(£) ® iifa) ® e fc (C)}, 

surface : S^ki^vX) = {k(0 ® ^(t?) ® e fe (C), e;(£) ® ® e fe (C), e<(0 ® e^) ® fc(Q}, 

volume : V$\(Z, V, C) = ® e^fa) ® e fe (C). 



Note that V^j. is indeed a 3-form, since ej(£) ® ej{rj) ® efc(C) = ei(0 £ j ( 7 ?) e fc(C) d£ A dry A d£. So 
the approximation spaces are spanned by combinations of Lagrange and edge basis functions, 

N,N,N 



A ,(Q;Co):=span{p^ fe } ' 

( m N,N,N ( , . ^N,N,N ( m > N,N,N 

AUQ; Cl ) := span (fe),)^ X { (4U} i=0 ,= life=0 >< ^(^L,^ 

r N.N,N f > N,N.N f N.N,N 

A^(Q;C 2 ) := s P an{(5S fe)l } x sp- { (^) 2 } i= w=1 x ^{(^.JsL^ 

A 3 ,(S;C 3 ):=span{^ 3 U , . ifc • 
t '-" j 1=1,1=1 ,fc=i 



Lagrange interpolation by itself does not guarantee a convergent approximation [32], but it 
requires a suitably chosen set of points, — 1 < £o < £i < • • • < £jv f_ L Here, the Gauss-Lobatto 
distribution is proposed, because of its superior convergence behaviour [35]. The convergence rates 
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of Lagrange and edge interpolants were obtained in [50] and are given by, 



(4.25) 
(4.26) 



\\a 



(0) 



^°>\\hao <C- 



pm- 1 

h 1 - 1 



- 7r fc &«|U» A i < C-^-I^^I^-iai, 



with I = min(p + l,m). For the variables vorticity, velocity and pressure in the VVP formulation 
of the Stokes problem, the /i-convergence rates of the interpolation errors become, 

\\u - n h uj\\ L2A n-2 = 0(h N+s ), \\u - 7r fc w|| H A-» = 0(h N ), 
(4.27) Hu-TThuH^-i =G(h N ), \\p--K h p\\ L * h « = 0{h N ), 

where s = 1 for n = 2 and s = for n > 2, and with N defined as in Section |4.3| Because of 
(3.20c I and (4.16), we have \\u — TTh.u\\ HA"- 1 = \\ u ~ 7I 7i m IIl 2 A"- 1 ■ 

4.4. Pointwise divergence-free discretization. One of the most interesting properties of the 
mimetic method presented in this paper, is that within our weak formulation, the divergence-free 
constraint is satisfied pointwise. This result follows from the three commuting properties with 
the exterior derivative, (4.10), (4.13) and (4.16), as was shown in LemmaJI] The corresponding 
commuting diagrams are repeated in the diagram below for the two dimensional case. 



K(Q;C ) — > A£(e ; Ci 



n 



V 



curl 



n 



K 



C°(D) 



-> C\D) 



-> C 2 {D) 



Note that by curl we refer to the two-dimensional variant, applied to a scalar, i.e. curlw = 
(dui/dy,—duj/dx) T , see also Example [lj and is also called the normal gradient operator, grad^, 
see [55] . 

In the following two examples we demonstrate the action of the exterior derivative on vorticity, 
S A°(Q re f; Co), and on the velocity flux, vtjp £ Al(Q le f; C\). Two dimensional reconstruction 
is based on tensor product construction of the one dimensional reconstruction function introduced 
above. 

Example 3 (Curl operator). Consider a flux z£ € hX(Q Te t;C\) with C\ outer-oriented, and 



where z 



(4.28) 



(i) 



dco 



(0) 



Then tc^ is expanded in the reference coordinates (£, n) as 

N N 



(0) 



i=0 j=0 



Apply the exterior derivative in the same way as in Lemma^ it gives 



dco 



(0) 



N N 

i=l j=0 
N N 



EE( 



N N 



i=0 j=l 



N N 



(4.29a) = - E E z lM0h(v) + E E iM^M, 

8=1 j=0 i=0 3=1 



where zf , 



and z?. 



-'i-lj 



can fee compactly written as z (!) = <W°\ rata 



w (o) g C°(D) and z^ 6 C X {D), or in matrix notation as z = E' 1,0 ^. TTiis relation is exact, 
coordinate free and invariant under transformations. 
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Example 4 (Divergence operator). Let £ A\(Q TC f ; C\) be the velocity flux defined as 

N N N N 

(4.30) 4 1 ' ^EE^'W+EE^W' 

■i=l j=o i=0 j=l 

Compare this to the velocity flux in ExamplelA p\9\ Then the change of mass, rnj^ £ A^(Q re f; C2 

(1) i— ' i— ' 

is equal to the exterior derivative of u h , 

N N 



(2) j (1) 



Ui 



(4.31) 



£5> 

t=i 3=1 

TV AT 

i=i 3=1 



1,3 -^j-Oe^Oej^)- 



where rriij = u^ j—u 



1,3+Wjj— Ui,j-i can 6e compactly written as = 6u^\ with u*- 1 ' € C 1 (Z)) 
and m^ 2 -* £ C 2 (D), or in matrix notation as m — E^'^u. Note that if the mass production 
is zero, as in our model problem (3.19b), the incompressibility constraint is already satisfied at 
discrete /cochain level. Interpolation then results in a pointwise divergence-free solution. 



5. Mixed formulation, boundary conditions and implementation 
Wc know how to discretize exactly the metric-free exterior derivative d (see Lemma [I] Sec- 



4.2 and the examples above), but it is less obvious how to treat the codiffcrcntial operator 
Fortunately, the two are directly related using L 2 -inner products as seen in (3.10). Therefore 



tion 
d* 

the deriva tion of t he mix ed formulation of the Stokes problem consists of two steps: 1). Multiply 



equations ( 3.20a )-( 3.20c I by the test functions r 



.(n-2) v (n 



using L 2 -inner products. 2). Use 



integration by parts, as in (3.10), to express the codifferentials in terms of the exterior derivatives 



and boundary integrals. The resulting mixed formulation of the Stokes problem becomes: 



Find (u>- 2 > 
{HA n - 2 x HA"- 1 



u (n-l) jp (n)) 



£ {HA n - 2 xHA n - 1 xL 2 A n }, given £ L 2 A n ~ 1 , for all (r 



(n-2) )W (n-l) g 



L 2 A n }, such that 



(5.1a) ( T (»- 2 ),o;("- 2 )) n - (dT^ n - 2 \u^) a = - / tr r^" 2 ) A tr W"" 1 ) 



(5.1b) (^"- 1 ),dc("- 2 )) n + (d t ;("- 1 ),p(")) n = (^"- 1 ),/("- 1 ))^+ / trv^Atr *p< n) , 



on 



(5.1c) 



( O ("),dn("- 1 )) fi =0. 



This mixed formulation is similar to those in [9j [3TJ [39] . The mixed formulation is well-posed, 
see [321 US]- The discrete problem is almost the same as the continuous problem, that is: find 



K 



£ A;;- 1 , for all (t<" 

xA" xA?}, such that ( 5.1a )-( 5.1c) hold. The discrete problem is also well-posed, because 



€ {Ar 2 xAj 



1 x AD, given /< n > 



-2) (n 



-2) (» 

{A£" 2 xA£ 

every subcomplex of a Hilbert complex is also a Hilbert complex, so if (HA, d) is a Hilbert complex, 
so is (A/j,d), and the projection operator from HA k (fl) to A^(fl;Ck) is bounded, see [S0J. A 
complete proof is given in . 

System (3.20) needs to be supplemented with boundary conditions on dtt. Their exists four 



possible types of boundary conditions as follows from the boundary integrals in the mixed formu- 
lation, (5.1 ). Subdivide the boundary into several parts, dfl = (L Tj, where Tj D Tj ; = for i 7^ j. 



Each part of the boundary can have one of the following four boundary conditions: 1. prescribed 
velocity (such as no-slip), 2. tangential velocity - pressure, 3. tangential vorticity - normal velocity, 
and 4. tangential vorticity - pressure boundary conditions. An overview is given in Table [l] 

From the implementation point of view we would like to mention that the L 2 inner products and 
boundary integrals are evaluated using Gauss-Lobatto quadrature, which is exact for polynomials 



Name 
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Type 



Normal velocity 
tangential velocity 


tr 


-l) ^ tr?/"" 1 ) 
tr * u y ' 


= 


u ■ n 


=> v ■ n = 
u ■ t 





essential 
natural 


Tangential velocity 
pressure 




tr * u v y 
tr 






U ■ t 
P 




natural 
natural 


Tangential vorticity 
normal velocity 


tr a/«- 
tr 


-2) ^ trr(™- 2 ) 
-l) ^ tru^"- 1 ) 


= 
= 


Q x t 
u ■ n 


=r* T X t — 

=>■ v ■ n = 


= 



essential 
essential 


Tangential vorticity 
pressure 


tr 


-2) ^ trr("- 2 ) 
tr *p( 2 ) 


= 


UJ X t 


r x r= 


= 


essential 
natural 



Table 1. Admissible boundary conditions for Stokes flow in vorticity- velocity- 
pressure formulation. 



up to order 27V — 1, [25]. The resulting system matrix is a saddle point system that is given by, 
(5.2) 



M (n-2) (' E («-l,n-2)j T M (r l -l) 








-Bi(*u) 






u 




M (n-l) f («-l) + B 2 (*p) 


M^Et"'™ -1 ) 




.p. 








The final system matrix is symmetric and only consists of L 2 inner product matrices for fc-forms, 
M( k > (also known as mass matrices), and incidence matrices, £i k ' k ~ x ) ^ that are directly obtained 
from the mesh topology, see p j!3| Coordinate transformations imposed by the pullback operator 
appear in the L inner products as a standard change of basis, see also |18j . The matrices Bi and 



B2 represent the boundary integrals in (5.1a) and (5.1b), and (*u) and (*p) are the tangential 
velocity and pressure boundary conditions imposed. A discussion on efficient solvers for symmetric 
indefinite systems that follow from saddle point problems can be found in [8] [72] . 



6. Numerical Results 

Now that all parts of the mixed mimetic method are treated, we can test the performance of the 
numerical scheme using a set of three test problems. The first one consists of an analytic solution 
on a unit square, where optimal /i-convergence and exponential p-convergence rates are shown for 
both Cartesian and curvilinear meshes for all combinations of boundary conditions. The second 
is a lid-driven cavity flow, where results are compared with a reference solution. Finally, Stokes 
flow around a cylinder moving with constant velocity in a channel is considered. 

6.1. Manufactured solution. The first test case addresses the convergence for h- and p-rcfincment 
of the mixed mimetic spectral element method applied to the Stokes model. The model problem 
is defined on the unit square f2 = [0, l] 2 , with Cartesian coordinates x := (x,y), with v — 1 and 
with the right hand side p 1 ' £ A 1 (ft) given by 



(i) 



(6.1a) 



dx + f x (x.) dy, 
— (-7T sin(-7ra;) cos(7ry) + 8tt 2 cos(27T2;) sin(27ry)) dx 
+ (7T cos(7re) sm(ny) — 8tt 2 sin(27rx) cos(27rj/)) dy. 
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This right hand side results in an exact solution for the vorticity w' ' € A°(f2), velocity flux 
g A 1 (fi), and pressure j^ 2 ) € A 2 (57) components of the Stokes problem, given by 

(6.1b) w (0) = w(x) = -47rsin(27ra;)sin(27r?/), 

i/ 1 ' = — f(x) da; + w(x) dy 
(6.1c) = — (cos(27rx) sin(27ry)) da; + (— sin(27ra;) cos(27ry)) dy, 

(6. Id) p*- 2 - 1 = p(x) dxAdy = (sin(7ra;) sin(7ry)) da;Ady. 

This testcase was discussed before in [38l [61] . Calculations were performed on both a Cartesian as 
well as a curvilinear mesh as shown in Figure 10 The map, (x,y) = $(£,77), used for the curved 
mesh is given by 

(6.2a) x(£, V) = h + Ht+I, sin «) sin(^)) , 

(6.2b) y% ,,) = | + §(,,+ | frinfrO sin(7T77)) . 




FIGURE 10. Examples of a Cartesian and a curvilinear mesh used in the conver- 
gence analysis. The meshes shown consist of 4 x 4 spectral elements, with for 
each element, N = 4. The element boundaries are indicated in red. 



Figure [TT] shows the /i-convergence and Figure 12 shows the p-convergence of the vorticity 
(4 0) G A° h {Q;C ), velocity e A l h {Q\ d) and pressure p% 5 € A 2 h (Q-C 2 )- For both figures, the 
results of the top row are obtained on Cartesian meshes and the results depicted underneath are 
obtained on curvilinear meshes. The errors for the vorticity and velocity are both measured in the 
L 2 A k - and f/A fe -norm, i.e. ||u;— aJh\\L 2 A°, w h||ifA°i an d \\u— m/i||l 2 a 1 , \\u— w/iIIjja 1 ) respectively. 
Because the divergence-free constraint is satisfied pointwise, the norm ||d(it — Wh)|U 2 A 2 is zero or 
machine precision, see Figure 13 and so the H A 1 -norm is equal to the L 2 A 1 -norm of the velocity, 
i.e., \\u — u/JIha 1 = \\ u ~ m /i||l 2 a 1 - This does not hold for the vorticity, since du/ ' € Al(Q;Ci) 
is again a function of sine and cosine functions. The norm ||d(w — ujf,)\\ L 2 A i converges one order 
slower than \\uj — 0Jh\\L 2 A - More details on the convergence behavior can be found in |49j . 



In Figure 11 the slope of the theoretical convergence rates, [49], are added which shows that 



/i-convergence rates are equal to the /i-convergence rates of the interpolation error (4.3|, on both 
Cartesian as well as curvilinear meshes. Figure [12] shows that exponential convergence rates are 
obtained on both types of meshes. 

It is important to remark is that these results are independent of the kind of boundary conditions 
used. This is shown in Table [2] This is an important result, because especially optimal conver- 
gence for the normal velocity - tangential velocity boundary condition is non-trivial in compatible 
methods, 4. The standard elements in compatible methods, the Raviart-Thomas elements, show 
only sub-optimal convergence for velocity boundary conditions, 0]. 
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Velocity 




Figure 11. Vorticity, velocity and pressure /i-convergence results of problem 



(6.1 1. Results in the top row correspond to Cartesian meshes, results in the 
bottom row are obtained on curvilinear meshes. All variables are tested on meshes 
with N = 2, 4, 6 and 8. 



Vorticity 




Velocity 





Figure 12. Vorticity, velocity and pressure p-convergence results of problem 



(6.1 1. Results in the top row correspond to Cartesian meshes, results in the 
bottom row are obtained on curvilinear meshes. All variables are tested on meshes 
with lxl, 2x2, 4x4 and 8x8 spectral elements. 



6.2. Lid-driven cavity Stokes. For many years, the lid-driven cavity flow was considered as one 
of the classical benchmark cases for the assessment of numerical methods and the verification of 
incompressible (Navier)-Stokes codes. The lid-driven cavity test case deals with a flow in a unit- 
square box with three solid boundaries and moving lid as the top boundary, moving with constant 
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A L -error 
▼ l) -error 
L -error 
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10 



10 

// 



10 



Figure 13. L 1 , L 2 and L°°-error of div u on the Cartesian mesh for discontinuous 
piecewise linear functions, N = 2. 



normal velocity 


tangential velocity 


vorticity 


vorticity 


convergence 


tangential velocity 


pressure 


normal velocity 


pressure 


rate 


4.0758e-01 


5.4293c-01 


5.4292e-01 


5.4292e-01 




1.9814e-01 


1.9738e-01 


1.9738e-01 


1.9738e-01 


1.46 


2.4893e-02 


2.4776e-02 


2.4776e-02 


2.4776e-02 


2.99 


3.1037e-03 


3.0954e-03 


3.0954e-03 


3.0954e-03 


3.00 


3.8738e-04 


3.8684e-04 


3.8684e-04 


3.8684e-04 


3.00 


4.8386e-05 


4.8352e-05 


4.8352e-05 


4.8351e-05 


3.00 


Table 2. This table shows the vorticity error — Uh\ 


i 2 A o obtained using the 



four types of boundary conditions described in Table [T] 
on a Cartesian mesh with N 
third order convergence. 



o aTir \ u — I I I -L J_ 

z, ctnu ii — 2' 4' 8' 16' 32 



The results are obtained 
, . All four cases show 



velocity equal to one to the right. Because of the discontinuities of the velocity in the two upper 
corners, the solution becomes singular at these corners, where both vorticity and pressure become 
infinite. Especially these singularities make the lid-driven cavity problem a challenging test case. 
For this test case a non-uniform 6x6 Cartesian spectral element mesh is used. Each spectral 



element consists of a Gauss-Lobatto mesh for N = 6, see Figure 14 The solutions of the vorticity, 



velocity, pressure and stream function are shown in Figure [14[ Also shown in Figure [14] is a plot of 
the divergence of velocity. It confirms a pointwise divergence-free solution up to machine precision. 
The results are in perfect agreement with those in |67j . 

Because in the mixed mimetic spectral element method no velocity unknowns are located at 
the upper corners - only velocity flux through edges is considered -, no special treatment is needed 
for the corner singularities, in contrast to many nodal finite-difference, finite-element and spectral 
element methods, [T71 [331 1HE1 EH] ■ This is due to the finite- volume like structure of the method, 
as explained in the section of algebraic topology. 

In Figure [15] the centerline velocities are plotted. Three different configurations are used, based 
on the same cell complex consisting of 9 x 9 2-cells: 

• left: 9x9 spectral elements with N = 1, resulting in piecewise constant approximations 
along the centerlincs, 

• middle: 3x3 spectral elements with N = 3, resulting in piecewise quadratic approxima- 
tions along the centerlines, 

• right: One global spectral element with N = 9, resulting in 8 th order polynomial approx- 
imations along the centerlines. 
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Wi)) 










FIGURE 14. Lid-driven cavity Stokes problem results. The top row from left to 
right shows the solution of the vorticity, velocity magnitude and pressure fields. 
The bottom row shows from left to right the solution of the stream function, the 
divergence of the velocity field and the 6 x 6, N = 6 mesh. 



Despite the low resolution, all approximations lay on top of those in [67] . 

Because of the tensor-product construction of discrete unknowns and basis-functions, an ex- 
tension to three dimensions is straightforward. A 3D lid-driven cavity is of interest because it 



not only contains corner singularities, but also line singularities. The left plot in Figure 16 shows 
slices of the magnitude of the velocity field in a three dimensional lid-driven cavity Stokes problem, 
obtained on a 2 x 2 x 2 element mesh with N — 8. The slices are taken at 10%, 50% and 90% 
of the y-axis. The right plot in Figure [16] shows slices of divergence of the velocity field. The 
solution at the symmetry plane coincides with the 2D results in Figure [Mj It confirms that also 
in three dimensions the mixed mimetic spectral element method leads to an accurate result with 
a divergence-free solution. 

The corner singularities can be made even more severe by sharpening the corners, as happens 
for a lid-driven cavity problem in a triangle. Figure [XT] shows the vorticity field and the velocity 
magnitude. On top of the velocity plot, stream function contours are plotted. The solutions are 
constructed on a 9 spectral element mesh with N — 9. A close-up of the stream function contours 



is shown in the rightmost plot in Figure 17 The stream function contours nicely show the first 
three Moffatt eddies (53) . 

6.3. Flow over a cylinder. The last test case considers the flow around a cylinder moving with 
constant velocity to the left, as defined in [26]. This testcase is mostly considered in the context 
of least-squares finite and spectral element methods, due to their moderate performance in case 
of large contraction regions, [2H [37] [55], mainly in terms of conservation of mass. 

The cylinder moves with unit velocity along the centerline of a narrow channel. The com- 



putational domain is defined as a rectangular box minus the cylinder, as shown in Figure 18 
Also visible in this figure are the 12 spectral elements in which the computational domain is di- 
vided. A transfinite mapping, [40] . is used to define the curved elements around the cylinder. 
Velocity boundary conditions of (u,v) = (1,0) are prescribed on the outer boundary and no-slip, 



JASPER. KREEFT AND MARC GERRITSMA 

9x9, N=1 3x3, N=3 1x1, N=9 




FIGURE 15. Horizontal (top) and vertical (bottom) centerline velocities are shown 
in blue for a very course mesh, 9x9 2-cells. From left to right the 9x9 2-cells 
are used in: 9x9 zeroth-order elements, 3x3 second-order elements and one 
eight-order element. In red the reference solution from [67]. 




FIGURE 16. Left: slices of magnitude of the velocity field of a three dimensional 
lid-driven cavity Stokes problem obtained on a 2 x 2 x 2 element mesh with N = 8. 
Right: slices of the divergence of velocity. Is confirms a divergence-free velocity 
field. 
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Figure 17. Lid-driven cavity Stokes flow in a triangle. Left the vorticity field, 
in the middle the velocity magnitude with stream function contours on top, and 
right the stream function contours of a close-up of the bottom corner, revealing 
the second and third Moffatt eddies. 



(u, v) = (0, 0), is prescribed along the boundary of the cylinder. Solution of the vorticity, velocity 



magnitude and pressure, together with streamlines are shown in Figure 18 




1.5 2 2.5 3 



FIGURE 18. Spectral element mesh (top left), magnitude of velocity (top right), 
vorticity (bottom left) and pressure (bottom right) for flow around a moving 
cylinder, on a 12 element, TV = 6 mesh. 



Next consider a control volume Q c consisting of the 6 elements in the domain —1.5 < x < 
0, 0.75 < y < 0.75. The control volume is chosen such that the ratio in size between inflow 
and outflow boundary is maximal. In this control volume conservation of mass should hold. 



Conservation of mass is expressed, by means of generalized Stokes theorem (3.6), in terms of a 
boundary integral as 



(6.3) 







du 



(1) £6 



.(!) 



From Section 4.4 and the results of the previous test cases we know that the solution of the velocity 
is divergence-free throughout the domain, independent of the chosen control volume. In Figure [19| 
a comparison is made for the horizontal velocity component u at the smallest cross-section above 
the cylinder, i.e. x — 0, 0.5 < y < 0.75, between the recently developed LSSCM, (35J, and our 
MMSEM method for N = 3,6, 12. Both methods use a similar mesh of 12 spectral elements. As 
can be seen from this figure, the MMSEM method performs already very well for N — 3, i.e. 
quadratic polynomial, where the LSSCM still fails for N = 6, i.e. sixth order polynomial. This is 
a direct consequence of the pointwise divergence-free discretization. 
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FIGURE 19. Horizontal velocity at smallest cross-section above the cylinder, on 
a 12 element mesh, for N — 3,6, 12. 

7. Conclusions and future aspects 

In this paper we presented the mixed mimetic spectral element method, applied to the vorticity- 
velocity-pressure formulation of Stokes model. At the heart lies the generalized Stokes theorem, 
which relates the boundary operator applied on an oriented geometric objects to the exterior deriv- 
ative, resembling the vector operators grad, curl and div, and the recently developed higher-order 
mimetic discretization for quadrilaterals and hexadrals, jSD]. The gradient, curl and divergence 
conforming method results in a point-wise divergence-free discretization of the Stokes problem, as 
was confirmed by a set of benchmark problems. These results also showed optimal convergence, 
independent of the type of boundary conditions on orthogonal and curved meshes. More on con- 
vergence behavior and error estimates is presented |49j . In the near future we plan to extend the 
method with structure-preserving /ip-refinement based on a compatible mortar element method. 
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